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DESCRIPTION OF PROGRAM 



PURPOSE 

There are many instances in which experimental data are available but 
no theoretical equation has been derived for expressing the relationship 
between two physical parameters. Sometimes the theoretical equation is 
known but is too complicated or cumbersome for use in mathematical evaluation. 
A polynomial is often used to approximate this relationship. If the coefficients 
of the polynomial are determined from the available data by the least squares 
procedure, the resulting equation can be used to evaluate the dependent 
variable with statistical validity from given values of the independent variable. 

This program generates an approximating polynomial by the least squares 
technique. The equation so derived contains as many terms as necessary to 
bring the standard error of the dependent variable within a range specified by 
the user. 

GENERAL DESCRIPTION 

Provision has been made for a maximum of 100 observations and a 
polynomial expansion up to the 15th degree. 

This program is based on the assumption that a set of experimental data 
can be fitted to a polynomial of the form: 

Y " AO + Alx +A2x 2 -+A15x 15 . 

The linear form Y s AO -+Alx is tried as a first approximation. The 
coefficients (A 1 a) are computed and the standard error of the dependent 
variable (Y) is compared to a predetermined maximum value. If the error 
is greater than the maximum value, the process is repeated, adding a term of 
the form AN xN (where N = 2, 3, . . . . , 15) until the error is within the 
tolerance. In each case the coefficients are computed by the least squares 
technique. The calculations utilize FORTRAN floating point arithmetic. A 
modified Gaussian elimination technique is used to solve the resulting set of 
linear equations. 

INPUT FORMAT 

All input to the program must conform to the specifications for fixed and 
floating point constants as outlined in the 1620 FORTRAN Manual. The input 
to the program consists of two types of records. The first record is read 
into the 1620 by Fortran Statement 1 which is "READ, N, TOL, LAST. " The 
next N records are the observations and are read into the 1620 by statement 
20 "READ, X(I), Y(I)" or if SW1 is on by statement 30 "READ, X(I), Y(I), 
W(I). " See page 1 of the appendix for explanation of the symbols used. The 
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exact format of these N records is dependent upon the setting of SW1. It is 
important to note that data are always entered with one observation to the 
record. 

If after obtaining a complete output, it is desired to enter a new TOL 
and/ or LAST, this can be accomplished through the console typewriter under 
control of SW4. See Operating Notes for details. 

OUTPUT 

All output is presently on the console typewriter. This can be easily 
modified by changing any of the PRINT statements to PUNCH statements. 
SW2 and SW3 control the format of the output. 

SW2 ON causes the following format to be printed for each order poly- 
nomial fitted: 

NORD TOL S2 N 

COEFFICIENT 

1 COEFFICIENT 



NORD COEFFICIENT 

SW3 ON causes the following additional output to be printed provided 
SW2 is also ON: 

X(l) Y(l) Y(l) Y(l) - Y(l) 

• • • • 

» • t • 

X(N) Y(N) Y(N) Y(N) - Y(N) 

All output will be obtained for the order polynomial being fitted if 
NORD = LAST or if S2 < TOL. If the output occurs under either or both of 
these conditions, the machine will halt after printing. See Operating Notes 
for switch settings before depressing START. 



MACHINE REQUIREMENTS 

Any IBM 1620 for which a FORTRAN compiler has been written. 
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THIS PROGRAM GENERATES AN APPROXIMATING POLYNOMIAL BY THE 
LEAST SQUARES TECHNIQUE. THE EQUATION SO DERIVED CONTAINS 
AS MANY TERMS AS NECESSARY TO BRING THE STANDARD ERROR OF 
THE DEPENDENT VARIABLE WITHIN A RANGE SPECIFIED BY THE 
USER, OR TO FIT A 15TH ORDER POLYNOMIAL. PRINTING OF 
INTERMEDIATE COEFFICIENTS AND THE PRINTING OF A TABULATION 
OF OBSERVED VS CALCULATED VALUES OF THE DEPENDENT VARIABLE 
ARE UNDER THE CONTROL OF PROGRAM SWITCHES AS IS THE 
INCLUSION OF WEIGHTING FACTORS. THE CALCULATIONS UTILIZE 
FLOATING ARITHMETIC WITH AN 8 DIGIT MANTISSA. IBM 1620, 
20K CORE, 1622 CARD READ-PUNCH. 



THIS PROGRAM AND ITS DOCUMENTATION WERE WRITTEN BY AN 
IBM EMPLOYEE. IT WAS DEVELOPED FOR A SPECIFIC PURPOSE AND 
SUBMITTED FOR GENERAL DISTRIBUTION TO INTERESTED PARTIES 
IN HOPE THAT IT MIGHT PROVE HELPFUL TO OTHER MEMBERS OF 
THE DATA PROCESSING COMMUNITY. THE PROGRAM AND ITS 
DOCUMENTATION ARE ESSENTIALLY IN THE AUTHORS ORIGINAL 
FORM. IBM SERVES AS THE DISTRIBUTION AGENCY IN 
SUPPLYING THIS PROGRAM. QUESTIONS CONCERNING THE USE OF 
THE PROGRAM SHOULD BE DIRECTED TO THE AUTHORS ATTENTION. 



Modifications or revisions to this program, as they occur, 
will be announced in the appropriate Catalog of Programs 
for IBM Data Processing Systems. When such an announce- 
ment occurs, users should order a complete new program 
from the Program Information Department. 



OPERATING NOTES 



OPERATING INSTRUCTIONS CARD SYSTEM 



1. Continue to Next Order Polynomial when Computer Halts after Printing. 

Turn SW4 ON 
Depress START 

When typewriter is activated, key in a value for TOL and LAST 

followed by a + , 
Turn SW4 OFF 
Depress RELEASE 
Depress START 

Note: If error is made while keying in this data, depress RELEASE and 

Start with SW4 ON. This will cause the typewriter to be activated 
and the whole record must be keyed in again. 

2. ERROR E7 OCCURS. 

You have probably divided by zero in the calculation of S2 between 
statement 160 and statement 161. The reason for this is that the 
order of the polynomial is within 1 of the number of observations. 
The coefficients calculated here will be correct, but the standard 
error S2 will be distorted. If you continue to the next order 
polynomial, ERROR F6 will occur as a result of taking the square 
root of a negative number when calculating S2. 
N 2 

3. If it is desired to compare V (Yi - Yi) rather than the present S2 

againBt TOL, delete the two FORTRAN statements immediately following 
statement 160 and re -assemble the program. 

4. Most questions about the program can be resolved by reference to the 
accompanying flow chart and FORTRAN listing. A complete analysis of 
the actual storage requirements of the program has not been made. The 
program did exceed memory by 1011 positions when assembled with the 
trace routine so there is very little space available for additions or lengthy 
modifications. 

5. A considerable amount of storage could be gained by deleting that portion 
of the program dealing with weight factors. 



I. Clear Memory 

RESET, INSERT, 310000300002, RELEASE, START 
INSTANT STOP 

II. Ready Typewriter for Normal Fortran Output 
Margins 10 and 95. 
Tabs every 17 positions. 

HI. Load Program 

Ready 1622 Read Hopper witti Program Deck 
RESET - 1620 
LOAD - 1622 

"LOAD DATA" will type out after loading is completed. 

IV. Ready Paper in Typewriter 

V. Place Data Cardtin 1622 Read Hopper 

VI. Set Console Switches 

OVERFLOW CHECK SWITCH - PROGRAM 
ALL OTHER CHECK SWITCHES - STOP 
PROGRAM SWITCHES: 



ON 

Weight Factors 

Print Polynomial 
Coefficients For Each 
Order Fitted 

Print Table of Observed vs 
Calculated After Listing 
Coefficients. 
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OFF 

No Weight Factors 

List Coefficients Only 
When NORD Z LAST or 
S2 ^ TOL. 

Print Table of Observed vs 
Calculated Only When 
NORD = LAST or S2 £ TOL. 

Read Next Set of Data 



Continue Present Problem 
(Normally off) 

* See pp. 3 for additional information on SW4. 



VII. Depress START ---- Or INSERT 4907500. 



COMPUTER 
TECHNOLOGY 



APPENDIX 
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SYMBOLS USED IN FORTRAN PROGRAM 



X(I) i^ 1 value of independent variable. 

Y(I) i^ 1 value of dependent variable. 

W(I) Weight factor of i** 1 observation. 

A(I, J) Matrix of coefficients of normal equations. 

N Number of observations. 

TOL Allowable value of standard error. 

LAST Highest order polynomial to be fitted. 
N 

Sum X(J) ^ (Wi) (Xi) J " 1 

Sum Y(J) jr (Yi) (Yi) (Xi) J " 1 
i= 1 

NORD Order polynomial fitted. 

SI Calculated value of Yi using Xi (also called Yi). 



S2 



1 



S3 Yi - *Yi 



N 

21 (Yi - YiV 

i Z 1 . 

(N - NORD - 1) 



POLYNOMIAL CURVE FITTING 
FLOW CHART 



Read, N, TOL, LAST 




OFF 



Read, X(I), 
*| Y(I) 

I = 1 to N 



Read, X(I), Y(I), W(I) 
I I 1 to N 



W(I) :1 
I: 1 to N 



Set up matrix for 
order 1. 




Accept, 
TOL, LAST 



OFF 



Solve for Polynomial 
Coefficients 



Calculate 
S2 




ON 



YES 



YES 





C POLYNOMIAL CURVE- F I ft I NG 

_C 

C .W R GRAVES IBM NEW ORLEANS ' 

C _ ^ _ 

D I MENS 10 H X ( 100),Y(100),A(ie,16),SUMX(31 ),SUMY(15),W(100) 

c 



1 READ, N.TOL, LAST 

DO hO 1=1 .N . 

IF(SENSE SWITCH 1)30,20 
20 READ.Xtl). Y(l) 

GO TO 40 
30 READ.Xfl ).Y(I ).W(I ) 



40 CONTINUE 

IFfSENSE SWITCH 1)70.50 



50 DO 60 l=t,N 

60_W(I)»1. 

'70 SUMX(1 )=0. 

SUMX(2)=0 1 __ 

SUMX(3)=0. 
SUMY(1)=0. 



SUMY(2)=0. 
DO 9 1=1 .N 



SUMX(1)=SUMX(1)+W(I ) 

S UMX (2) =SUMX(2)+W( I )*X( I ) 

SUMX(3 )=SUMX(3 )+W( I )*X( I )*X( 1 
SUMY( 1 )=SUMY(1 )+W( I )*Y( I ) 
90 SUMY(2)=SUMY(2)+W(I )*X( I )*Y(I ) 



N0RD=1 



91 I F( SENSE SWITCH 4)92,93 

92 ACCEP T. TOL. LAST 



93 L=N0RD+1 
KK=L+1 



DO 101 1=1 ,L 
DO 100 J=1 f L 



I K=J-1+I 
100 A( I .J)-SUMX(jK) 



101 A(i,KK)=SUMY(l ) 
DO 140 1=1 .L 



A(KK, I )=-1 . 
KKK=I+1 



DO 110 J=KKK,KK \^ 



110 A(KK f J)=0. 



C-1./A(1,l) 

DO 120 I 1=2. KK 



DO 120 J=KKK,KK 
120 API , J )=A( I I .J)-A(1.J)*A(U.I)*C 
DO 140 I 1 = 1 ,L 
DO 140 J=KKK.KK 



140 A(l l,J)=A(l 1+1, J) 

S2=0. 

DO 160 J=1 ,N 

S1_=jO._ 



S1=S1+A(1 ,KK) 
DO 150 1=1 ,NORD 



150 S1=S1+A(I+1,KK)*X(J)**I 
160 S2=S2+(S1-Y(J))*(S1-Y(J)) 



3=THT 

S2=(S 2VB)**.5 



I FTSENSE SWITCH 2)163,161 
1 61 IF(N0RD-LAST)162.163.162 



162 i F(S2-T0L)163 ,163 ,171 

163 PRINT, " n \, 
PRINT, NORD, TOL, S2, N 



DO 164 1=1, L 
J=l-1 

164 PRINT, J ,A( I ,KK) ~ ~ 

. . ...I F( SENSE SWITCH _J 11 62*165 

165 I F ( NORD-LAST) 1 667167 ,t66 " ~ 

166 I Ff S2-TQL) 167 ' 167 171 

167 DO 169 1 = 1, N ~ ' 

S1=0. 

S1=A(1,KKT " — 

D0_168„J=1 ,NORD 

168 Si=S1+A"( J+1 ,KKT*X( I )**J : 

S3=Yfl)-S1 

169 PRINT, X(l),Y(l ),S1,S3 " 

LF(NQRD- lAST) 170. 173.1 73 

170 IF(S2-TOL)173,173,171 : " 

_LZ1_N0RD=N0 RD+1 - 

J=2*N0RD " ' " ~ 

suMx(j)=n. 

SUMX(J+1)=0. ~ " — 

SUMY( NQRD+li =0. 

DO 172 1=1 ,N ~ ~ — 

S UMXf J)=SUMX(J)+X( I )**(J-1 )*wf | ) 

SUMXU+1 )=SUMX(J+1 )+X(l )**J*W( I) ~ 

172 GQ M TQ N ° RD " t ' 1 )=SUMY(N0 ^ +1 )+Y ( 1 )*XC I )**NORD*Wf I ) 

173 PAUSE " 

IF(SENSE SWITCH 4)171,1 ~ 

END 



0- 
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SAMPLE PROBLEM 1 



NO WEIGHT FACTORS 



n 



-5.0000000 
=4^0000000- 
-3.0000000 
^0000000- 
-1.0000000 
000000. 
1 .0000000 
2.0000000 



*i, . 3.0000000 
i ' 4.oooooon 



5.0000000 



_,1..0000000E=03_ 
2.0000000 

-^-99999999 

7.0000000 
_fi-J)000000 

5.0000000 
_AJ)000000 

3.0000000 
__L_Q00000Q . 



1 .0000000 
_JIQ00Q0.QD_ 
-1 .0000000 

-?.nonnnno 



-4»7l4o470E=a8- 



7.0000000 

.6.. ooooooo- 

5.0000000 

Juooooooo_ 

2.9999999 
JUOOOOOOCL 

1 .0000001 
_J&Q0D000CL 



-1 .0000000 

-?.oooooon 



41 



.00000000 
^00000000- 

.00000000 
-J)0000000_ 



1 .0000000E-07 
_jQ000000Q 



-3.0000000 



-3.0000000 



-1 .0000000E-07 
.onnnonn n 
.00000000 

.nnnnonno 



.00000000 



c 

c 
c 
c 
o 

o 



SAMPLE PR OBL EM 2 NO WEIGHT FACTORS 



00000000 .137096JA 35- 

.67653961* 
k»5094777E-Q2 



.00000000 6.8247925E-02 35 
-88567484 v - 



-5.8891750E-O2 
7.3 669823 E=Q3_ 



.00000000 2.2010479E-02 35 

U01 01370__ 

-.16878278 

2^7903 1 »53Ej = 02_ 



-9. 6864043 E -04 

. 00.000000 1.266765 4E - Q2 35 

1.0527606 
--Z2962000. 



'4.8020588E-02 
-3lJ£7064_QE=Q3_ 



7.51 53401 E-05 

.Q0BP0000 1. 252851 8E-02 35 



1.045601' 

-^21550438 

4T1142057E-02 
-1 .893 9045 E -03 



-2.2981 124T^05 
2 .658 6468E-06 



00000000" 1 .4404736E-02 3T 

1 .0549349 

-.23916686 

5_.6848405E r 02 

-6.1416738E-03 

5.1292009E-0 4 



-27883641 6E -05 

6 .967663 8E-07 ' 

.00000000' 1 .166662 1E-02 3T 
1 .041958 2 



-.20453217 
3. 3409581 E-02 



3.1708200E-05 

-U81 865 13 E^Q4 

-1 .2040700E-OS 

1 . 074092 7 E -06 

-3 .61 92355T=08 
. 00000000 1 .2140386E-02 35 



1 .0428216 

-.20598632 

3 .3852269-^2^ 

6.47968l3E-0<j ' 

-27 06 64440E -64 

7.8117300E-0 7 

1 .1782136E-06 "~. 

-5 t 4l57920E-08 PARTIAL OUTPUT 

575021 779£ r T5 

. 0000000 . ■ 1.3673123E-02 35 

1Td537794 : ~ : 

-.22739792 ; 

4.3658105E-02 
-1.0574389E-03 
-33231047E-04' 



3.7879753E-0 5 IK 



47B058560E-07 
-3.42033 67E-07 

-5.2506785E-10 



-=5^0000000- 
-4.0000000 
^=3.. 0000000. 
-2.0000000 
-1 .0000000- 
.00000000 
1 .0000000- 
2.0000000 
s.oononnn 
if. 0000000 
5.0000000 



^30000000 
3.0090909 

-J. 8918182 
.30000000 

-U0160840 
1.8918182 

—19930069 
.30000000 

_U0 160840 
1 .0097709 
J3930069 



4.9553224E-02 

^.2500000 

-3.0000000 

.si .55 00000 

-.60000000 

.15000000 

1 .0000000 
2.7500000 



4.3000000 
_LJ5J>Q000_ 



1 t .400000 
17.250000 



MP LB PROBLEM 3 NO WEIGHT FACTORS 



-2-^3943-79- 



4X 



1 .3773300 



11 



3.1790805E-02 



11 



^244406 2_ 
-3.0055949 
^:U5574595_ 
-.60268089, 
__a5606057— 
1.0160840 
2.27-4708Z_ 
4.2292542 
-7.1ZZ0399— 



11.415384 
17.241608 



^=5_5938000E^03_ 
5.5949000E-03 

_7-4595000E-03_ 
2.680890OE-O3 

^=6^O6O57OOE^03- 

-1 .6084000E-02 



7.0745800E-02 
j=Z±1Q339QQZ=QZ- 



-1 .5384000E-02 
„ 8 . 39 200 0E-Q 3 



